Cream of the Crop 20/Cream of the Crop 20 (Terry Blount) (1996).iso/os2/xdsn217.zip/SAMPLES/SIMPLE/linnew.mod

Text File | 1996-07-10 | 26KB | 903 lines

<* IF __GEN_X86__ THEN *>
  <*+NOPTRALIAS*>
  <*-SPACE*>
  <*-GENHISTORY*>
  <*+DOREORDER*>
<* END *>
<* ALIGNMENT="4"*>
<*+PROCINLINE*>
<*-CHECKINDEX*>
<*-CHECKRANGE*>
<*-CHECKNIL*>
<*-IOVERFLOW*>
<*-COVERFLOW*>
<*-GENDEBUG*>
<*-LINENO*>

<* IF NOT DEFINED (DP) THEN *>
<* NEW DP+ *>
<* END *>

MODULE linnew;

IMPORT SYSTEM;
IMPORT SysClock;
FROM STextIO IMPORT WriteString, ReadString, WriteLn, SkipLine;
FROM <* IF DP THEN *> SLongIO <* ELSE *> SRealIO <* END *> IMPORT WriteFixed;
FROM SWholeIO IMPORT WriteInt;
FROM WholeStr IMPORT StrToInt;
FROM ConvTypes IMPORT ConvResults;
FROM Storage IMPORT ALLOCATE, DEALLOCATE;

(*
**
** LINPACK.C        Linpack benchmark, calculates FLOPS.
**                  (FLoating Point Operations Per Second)
**
** Translated to C by Bonnie Toy 5/88
**
** Modified by Will Menninger, 10/93, with these features:
**  (modified on 2/25/94  to fix a problem with daxpy for
**   unequal increments or equal increments not equal to 1.
**     Jack Dongarra)
**
** - Defaults to double precision.
** - Averages ROLLed and UNROLLed performance.
** - User selectable array sizes.
** - Automatically does enough repetitions to take at least 40 CPU seconds.
** - Prints machine precision.
** - ANSI prototyping.
**
** To compile:  cc -O -o linpack linpack.c -lm
**
**
*)

<* IF DP THEN *>

CONST PREC      = "Double";
      BASE10DIG = 14; (* DBL_DIG; *)

TYPE FLOAT = LONGREAL;

<* ELSE *>

CONST PREC      = "Single";
      BASE10DIG = 7; (* FLT_DIG; *)

TYPE FLOAT = REAL;

<* END *>

TYPE
    FLOATARRAY    = ARRAY [0..10000000] OF FLOAT;
    INTARRAY      = ARRAY [0..10000000] OF INTEGER;
    FLOATARRAYPTR = POINTER TO FLOATARRAY;
    INTARRAYPTR   = POINTER TO INTARRAY;

(*
** Constant times a vector plus a vector.
** Jack Dongarra, linpack, 3/11/78.
** ROLLED version
*)

PROCEDURE daxpy_r (n : INTEGER; da : FLOAT;
                   VAR dx : FLOATARRAY; incx : INTEGER;
                   VAR dy : FLOATARRAY; incy : INTEGER);

VAR i,ix,iy : INTEGER;
BEGIN
    IF n <= 0   THEN RETURN END;
    IF da = 0.0 THEN RETURN END;

    IF (incx <> 1) OR (incy <> 1) THEN

        (* code for unequal increments or equal increments <> 1 *)

        ix := 1;
        iy := 1;
        IF incx < 0 THEN ix := (-n+1)*incx + 1; END;
        IF incy < 0 THEN iy := (-n+1)*incy + 1; END;
        FOR i := 0 TO n-1 DO
            dy[iy] := dy[iy] + da*dx[ix];
            INC (ix, incx);
            INC (iy, incy);
        END;
        RETURN;
    END;

    (* code for both increments equal to 1 *)

    FOR i := 0 TO n-1 DO
        dy[i] := dy[i] + da*dx[i];
    END;
END daxpy_r;

(*
** Forms the dot product of two vectors.
** Jack Dongarra, linpack, 3/11/78.
** ROLLED version
*)

PROCEDURE ddot_r (n : INTEGER; VAR dx : FLOATARRAY; incx : INTEGER;
                               VAR dy : FLOATARRAY; incy : INTEGER) : FLOAT;

VAR dtemp : FLOAT;
    i,ix,iy : INTEGER;
BEGIN
    dtemp := 0.0;

    IF n <= 0 THEN RETURN 0.0; END;
    IF (incx <> 1) OR (incy <> 1) THEN

        (* code for unequal increments or equal increments <> 1 *)

        ix := 0;
        iy := 0;
        IF incx < 0 THEN ix := (-n+1)*incx; END;
        IF incy < 0 THEN iy := (-n+1)*incy; END;
        FOR i := 0 TO n-1 DO
            dtemp := dtemp + dx[ix]*dy[iy];
            INC (ix, incx);
            INC (iy, incy);
        END;
        RETURN dtemp;
    END;

    (* code for both increments equal to 1 *)

    FOR i:=0 TO n-1 DO
        dtemp := dtemp + dx[i]*dy[i];
    END;
    RETURN dtemp;
END ddot_r;


(*
** Scales a vector by a constant.
** Jack Dongarra, linpack, 3/11/78.
** ROLLED version
*)

PROCEDURE dscal_r (n : INTEGER; da : FLOAT;
                   VAR dx : FLOATARRAY; incx : INTEGER);

VAR i : INTEGER;
BEGIN
    IF n <= 0 THEN RETURN END;
    IF incx <> 1 THEN

        (* code for increment not equal to 1 *)

        FOR i := 0 TO n-1 DO
            dx[i*incx] := da*dx[i*incx];
        END;
        RETURN;
    END;

    (* code for increment equal to 1 *)

    FOR i := 0 TO n-1 DO
        dx[i] := da*dx[i];
    END;
END dscal_r;

(*
** constant times a vector plus a vector.
** Jack Dongarra, linpack, 3/11/78.
** UNROLLED version
*)

PROCEDURE daxpy_ur (n : INTEGER; da : FLOAT;
                    VAR dx : FLOATARRAY; incx : INTEGER;
                    VAR dy : FLOATARRAY; incy : INTEGER);

VAR i,ix,iy,m : INTEGER;
BEGIN
    IF n <= 0   THEN RETURN END;
    IF da = 0.0 THEN RETURN END;

    IF (incx <> 1) OR (incy <> 1) THEN

        (* code for unequal increments or equal increments <> 1 *)

        ix := 1;
        iy := 1;
        IF incx < 0 THEN ix := (-n+1)*incx + 1; END;
        IF incy < 0 THEN iy := (-n+1)*incy + 1; END;
        FOR i := 0 TO n-1 DO
            dy[iy] := dy[iy] + da*dx[ix];
            INC (ix, incx);
            INC (iy, incy);
        END;
        RETURN;
    END;

    (* code for both increments equal to 1 *)

    m := n MOD 4;
    IF m <> 0 THEN
        FOR i := 0 TO m-1 DO
            dy[i] := dy[i] + da*dx[i];
        END;
        IF n < 4 THEN
            RETURN;
        END;
    END;
    FOR i := m TO n-1 BY 4 DO
        dy[i]   := dy[i] + da*dx[i];
        dy[i+1] := dy[i+1] + da*dx[i+1];
        dy[i+2] := dy[i+2] + da*dx[i+2];
        dy[i+3] := dy[i+3] + da*dx[i+3];
    END;
END daxpy_ur;

(*
** Forms the dot product of two vectors.
** Jack Dongarra, linpack, 3/11/78.
** UNROLLED version
*)

PROCEDURE ddot_ur (n : INTEGER;
                   VAR dx : FLOATARRAY; incx : INTEGER;
                   VAR dy : FLOATARRAY; incy : INTEGER) : FLOAT;

VAR dtemp : FLOAT;
    i,ix,iy,m : INTEGER;
BEGIN
    dtemp := 0.0;

    IF n <= 0 THEN RETURN 0.0; END;

    IF (incx <> 1) OR (incy <> 1) THEN

        (* code for unequal increments or equal increments != 1 *)

        ix := 0;
        iy := 0;
        IF incx < 0 THEN ix := (-n+1)*incx; END;
        IF incy < 0 THEN iy := (-n+1)*incy; END;
        FOR i := 0 TO n-1 DO
            dtemp := dtemp + dx[ix]*dy[iy];
            INC (ix, incx);
            INC (iy, incy);
        END;
        RETURN dtemp;
    END;

    (* code for both increments equal to 1 *)

    m := n MOD 5;
    IF m <> 0 THEN
        FOR i := 0 TO m-1 DO
            dtemp := dtemp + dx[i]*dy[i];
        END;
        IF n < 5 THEN
            RETURN dtemp;
        END;
    END;
    FOR i := m TO n-1 BY 5 DO
        dtemp := dtemp + dx[i]*dy[i] +
                         dx[i+1]*dy[i+1] +
                         dx[i+2]*dy[i+2] +
                         dx[i+3]*dy[i+3] +
                         dx[i+4]*dy[i+4];
    END;
    RETURN dtemp;
END ddot_ur;


(*
** Scales a vector by a constant.
** Jack Dongarra, linpack, 3/11/78.
** UNROLLED version
*)

PROCEDURE dscal_ur (n : INTEGER; da : FLOAT;
                    VAR dx : FLOATARRAY; incx : INTEGER);

VAR i, m : INTEGER;
BEGIN
    IF n <= 0 THEN RETURN END;
    IF incx <> 1 THEN

        (* code for increment not equal to 1 *)

        FOR i := 0 TO n-1 DO
            dx[i*incx] := da*dx[i*incx];
        END;
        RETURN;
    END;

    (* code for increment equal to 1 *)

    m := n MOD 5;
    IF m <> 0 THEN
        FOR i := 0 TO m-1 DO
            dx[i] := da*dx[i];
        END;
        IF n < 5 THEN RETURN END;
    END;
    FOR i := m TO n-1 BY 5 DO
        dx[i]   := da*dx[i];
        dx[i+1] := da*dx[i+1];
        dx[i+2] := da*dx[i+2];
        dx[i+3] := da*dx[i+3];
        dx[i+4] := da*dx[i+4];
    END;
END dscal_ur;

(*
** Finds the index of element having max. absolute value.
** Jack Dongarra, linpack, 3/11/78.
*)

PROCEDURE idamax (n : INTEGER; VAR dx : FLOATARRAY; incx : INTEGER) : INTEGER;

VAR dmax         : FLOAT;
    i, ix, itemp : INTEGER;

BEGIN
    IF n < 1 THEN RETURN -1 END;
    IF n = 1 THEN RETURN  0 END;
    IF incx <> 1 THEN

        (* code for increment not equal to 1 *)

        ix := 1;
        dmax := A